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ABSTRACT 

Using three-dimensional magnetohydrodynamics simulations, we investigate general properties of a blast 
wave shock interacting with interstellar clouds. The pre- shock cloudy medium is generated as a natural conse- 
quence of the thermal instability that simulates realistic clumpy interstellar clouds and their diffuse surrounding. 
The shock wave that sweeps the cloudy medium generates a turbulent shell through the vorticity generations 
that are induced by shock-cloud interactions. In the turbulent shell, the magnetic field is amplified as a result 
of turbulent dynamo action. The energy density of the amplified magnetic field can locally grow comparable to 
the thermal energy density, particularly at the transition layers between clouds and the diffuse surrounding. In 
the case of a young supernova remnant (SNR) with a shock velocity ^10^ km s~\ the corresponding strength 
of the magnetic field is approximately 1 mG. The propagation speed of the shock wave is significantly stalled 
in the clouds because of the high density, while the shock maintains a high velocity in the diffuse surrounding. 
In addition, when the shock wave hits the clouds, reflection shock waves are generated that propagate back 
into the shocked shell. From these simulation results, many observational characteristics of a young SNR RX 
J17 13.7-3946 that is suggested to be interacting with molecular clouds, can be explained as follows: The re- 
flection shocks can accelerate particles in the turbulent downstream region where the magnetic field strength 
reaches ImG, which causes short-time variability of synchrotron X-rays. Since the shock velocity is stalled 
locally in the clouds, the temperature in the shocked cloud is suppressed far below 1 keV. Thus, thermal X- 
ray line emission would be faint even if the SNR is interacting with molecular clouds. We also find that the 
photon index of the Tr^-decay gamma rays generated by cosmic-ray protons can be 1.5 (corresponding energy 
flux is vFi, oc z/^^), because the penetration depth of high-energy particles into the clumpy clouds depends on 
their energy. This suggests that, if we rely only on the spectral study, the hadronic gamma-ray emission is 
indistinguishable from the leptonic inverse Compton emission. We propose that the spatial correlation of the 
gamma-ray, X-ray, and CO line emission regions can be conclusively used to understand the origin of gamma 
rays from RX J1713.7-3946. 

Subject headings: acceleration of particles — gamma rays: ISM — ISM: supernova remnants — magnetic 
fields — shock waves — turbulence — X-rays: individual (RX J1713.7-3946) 



1. INTRODUCTION 

Supernova remnants (SNRs) are believed to be the sites of 
Galactic cosmic-ray acceleration through a diffusive shock 
acceleration mechanism (DSA; Bell 1978, Blandford & Os- 
triker 1978, Blandford & Eichler 1987), and multi- wavelength 
nonthermal emissions from SNRs caused by accelerated par- 
ticles have been detected (e.g., Koyama et al. 1995, Aharo- 
nian et al. 2008). However the detailed process of DSA and 
the emission mechanism of SNRs are still a matter of debate. 
Recent observations have emphasized the importance of the 
interaction between SNRs and interstellar clouds. The Fermi 
Gamma-ray Space Telescope revealed gamma-ray emissions 
from middle-aged SNRs interacting with molecular clouds 
(Abdo et al. 2009, 2010a, 2010b, 2010c). In addition, in a 
young SNR, RX J1713.7-3946 (or G347.3-0.5), spatial cor- 
relation between molecular clouds and X-/gamma-ray emis- 
sions has been reported (Fukui et al. 2003, 2008, 2011, 
Moriguchi et al. 2005, Sano et al. 2010). 

When we consider the interaction between a SNR and in- 
terstellar clouds, we should take into account the highly inho- 
mogeneous structure of clouds (see also Laming 2001a, b for 
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the effects of circumstellar density inhomogeneity). Recent 
numerical simulations have shown that interstellar clouds are 
formed as a complex of clumps fragmented by the thermal 
instability that are embedded in diffuse gas (Koyama & Inut- 
suka 2002; Hennebelle et al. 2008; Inoue & Inutsuka 2008; 
2009; Banerjee et al. 2009; Heitsch et al. 2009; Vazquez- 
Semadeni et al. 2006, Audit & Hennebelle 2010). Further- 
more, in molecular clouds, supersonic turbulence is always 
observed as a supra-thermal line width of molecular line emis- 
sions (e.g., Larson 1981, Heyer & Brunt 2004) that inevitably 
generate highly inhomogeneous structures by shock compres- 
sions (see, e.g., MacLow & Klessen 2004). By using two- 
dimensional magnetohydrodynamics (MHD) simulations, In- 
oue, Yamazaki & Inutsuka (2009) demonstrated that the in- 
teraction between a strong shock wave and a cloudy inhomo- 
geneous medium generates a turbulent SNR shell. The turbu- 
lence induced by the shock-cloud interactions amplifies the 
magnetic field through turbulent dynamo actions (see also, 
Balsara et al. 2001; Giacalone & Jokipii 2007), which po- 
tentially account for the year- scale short- time variability of 
X-rays in RX J1713.7-3946 (Uchiyama et al. 2007) as well 
as the spatial scale of the regions with the short-term variation. 
In the follow-up study by Inoue, Yamazaki & Inutsuka (2010), 
a similar but more realistic three-dimensional simulation was 
performed in which it was found that the shock-cloud interac- 
tions cause reflected shock waves in the SNR that were shown 
to reproduce the broken power-law cosmic-ray spectrum ob- 
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served in the above-mentioned middle-aged SNRs (see also 
Malkov et al. 2010, Ohira et al. 2010, Uchiyama et al. 2010). 
The shock velocity examined in Inoue, Yamazaki & Inutsuka 
(2010) is much smaller than that expected in young SNRs be- 
cause we discussed middle-aged SNRs and focused only on 
the physical properties of the reflected shock waves. In this 
paper, we study the shock-cloud interaction in more detail 
and examine its influences especially on young SNRs aged 
- 1000 years such as RX J1713.7-3946. 

A number of observational and theoretical efforts have 
been devoted to understand the emission mechanisms of RX 
J1713.7-3946 (Enomoto et al. 2002, Butt et al. 2002, Aha- 
ronian et al. 2006, 2007, Uchiyama et al. 2007, Tanaka et 
al. 2008, Kats & Waxman 2008, Plaga 2008, Berezhko & 
Volk 2008, Acero et al. 2009, Fang et al. 2009, Morlino 
et al. 2009, Yamazaki et al. 2009, ZirakashviH & Aharo- 
nian 2010, ElHson et al. 2010, Zhang & Yang 2011, Abdo 
et al. 2011). Based on the recent gamma-ray observation 
using the Fermi space telescope, it is suggested that the spec- 
tral energy distribution prefers the emission mechanism of the 
inverse Compton scattering of the cosmic microwave back- 
ground (CMB) photons by high-energy electrons accelerated 
by the shock wave (Abdo et al. 2011). Theoretical mod- 
eling of RX J17 13.7-3946 also supports the leptonic origin 
of the gamma-ray emission; Ellison et al. (2010) used one- 
dimensional hydrodynamics simulations of a supernova blast 
wave to demonstrate that the hadronic model of gamma-ray 
emission fails to reproduce observed X-ray emission due to 
an overproduction of thermal X-ray line emission. How- 
ever, the leptonic model of gamma-ray emission from RX 
J 17 13.7-3946 also has problems. If the gamma-ray emis- 
sion is leptonic, the magnetic field strength should he B ^ 10 
/iG from the ratio of synchrotron X-ray and inverse Compton 
gamma-ray fluxes, while the short-time variability observed 
in some X-ray bright regions indicates that the magnetic field 
strength can be as large as 5 ~ 1 mG, and the thickness of the 
X-ray filaments indicates B ^ 100 jiG around the shock front 
(Ballet 2006, Uchiyama et al. 2007, Tanaka et al. 2008, Acero 
et al. 2009), which is similar to other young SNRs (Vink & 
Laming 2003, Bamba et al. 2003, 2005). 

The lack of thermal X-ray line emission from RX 
J 17 13.7-3946 is apparently very crucial when we develop a 
shock-cloud interaction model of young SNRs, because the 
shocked dense clouds are thought to emit numerous thermal 
X-ray lines. In this paper, however, using three-dimensional 
MHD simulations, we show that the clouds shocked by a su- 
pernova blast wave do not emit thermal X-ray lines, because 
the transmitted shock wave stalls in the dense clouds. This 
may resolve the problem that the thermal X-ray line emis- 
sion is substantially suppressed in RX J17 13.7-3946 despite 
the suggested interaction with dense molecular clouds. We 
also show that the photon index of the hadronic gamma rays 
can be p-l/2= 1.5 for p = 2, where p is the spectrum in- 
dex of accelerated protons, which is consistent with the recent 
gamma-ray observation by Abdo et al. (201 1), if the interact- 
ing molecular clouds are clumpy as demanded from theoreti- 
cal arguments, numerical simulations, and observational facts 
mentioned above. This indicates that it is difficult to distin- 
guish the leptonic and hadronic gamma-ray emissions from 
the spectral study alone. We propose the spatial correlation 
of the gamma-ray. X-ray, and CO line emission regions as a 
conclusive tool to understand the origin of gamma rays from 
RX J1713.7-3946. 

The organization of the paper is as follows: In §2, we briefly 
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Fig. 1 . — Thermal equilibrium curve in typical ISM in which cooling and 
heating are in balance (thick solid). A fitting cooling/heating function given in 
Koyama & Inutsuka (2002) is used to calculate the curve. The dotted region 
is the thermally unstable region in which the Balbus instability criterion eq. 
{l] is satisfied. Dashed lines are isotherms ofT= 10^, 10^, and 10^ K. 

explain why we consider inhomogeneous clouds and their sur- 
rounding and provide numerical settings of our simulations. 
The results of the simulations are shown in §3. In §4, we dis- 
cuss implications of the simulations and their application to 
the young SNR RX J1713.7-3946. Finally in §5, we summa- 
rize our findings. 

2. SETUP OF SIMULATIONS 

2.1. Thermal Instability in the ISM 

In this paper, we consider inhomogeneous interstellar 
clouds as a pre-shock medium. In the following, we briefly 
explain how inhomogeneity is imprinted in interstellar clouds. 
It is widely known that the interstellar medium (ISM) is an 
energy-open system due to the radiative cooling and heating 
that make the ISM a thermally bistable medium (Field et al. 
1969; Wolfire et al. 1995; 2003; Koyama & Inutsuka 2000). 
Fig. T] shows the thermal equilibrium curve in a typical ISM 
in which cooling and heating are in balance (Koyama & In- 
utsuka 2002). Two roughly isothermal states with T ^ 10^ 
K and 10"^ K correspond, respectively, to the phases of in- 
terstellar cloud and diffuse intercloud gas that are connected 
by a thermally unstable equilibrium. Below the equilibrium 
curve, cooling dominates heating, and above it, heating dom- 
inates cooling. The thermal instability is the most promising 
formation mechanism of interstellar clouds that is driven by 
runaway cooling. The instability criterion is given by the con- 
dition (Balbus 1995): 



d_ 
df 



<0, 



(1) 



where jC(n^ T) is the net cooling rate per unit mass. The dotted 
region in Fig. [T]is the thermally unstable region in which den- 
sity inhomogeneities grow exponentially toward stable equi- 
libria (see. Field 1965, Schwarz et al. 1972, and Koyama & 
Inutsuka 2000 for linear stability analyses under various dy- 
namical conditions). From Fig. [T] it is clear that, during any 
formation process of clouds from the diffuse intercloud gas, 
the gas always experiences thermal instability. In a typical 
ISM, the timescale of the thermal instability is given by the 
cooling timescale that is ^cooi ~ 1 Myr, and the most unsta- 
ble scale of the thermal instability is /ti ~ 1 pc. The non- 
linear growth of the thermal instability toward the cold phase 
(condensation) generates a cold clump whose scale is much 
smaller than /ti. This is the reason why recent numerical sim- 
ulations of molecular and HI clouds formation show the gen- 
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eration of clouds as a complex of cloudlets (Koyama & Inut- 
suka 2002; Hennebelle et al. 2008; Inoue & Inutsuka 2008; 
2009; Banerjee et al. 2009; Heitsch et al. 2009; Audit & Hen- 
nebelle 2010). It is also known from the comparisons of the 
growth timescales of various hydrodynamic instabilities that 
the thermal instability dominates the dynamics of cloud for- 
mation (Heitsch et al. 2008). The simulations of the global 
ISM under the influence of supernovae have also pointed out 
that the multiphase structure in the ISM is regulated by the 
dynamical process of the thermal instability (MacLow et al. 
2005, de Avillez & Breitschwerdt 2005). Observationally, 
Sakamoto & Sunada (2003) found, at the envelope in the Tau- 
rus molecular cloud, that the cloud is indeed composed of the 
small-scale clumps where the thermal instability plays a role. 

In molecular clouds, there is an additional generator of den- 
sity inhomogeneity other than the thermal instability. Molec- 
ular line emissions from clouds are always observed with 
supra-thermal line widths that are considered to be the out- 
come of supersonic turbulence (Larson 1981, Heyer & Brunt 
2004, Mac Low & Klessen 2004). The supersonic turbulence 
inevitably introduces density fluctuations due to the shock 
compressions of converging flows and rarefactions by diverg- 
ing flows, even if we artificially set up a uniform molecular 
cloud initially (see, e.g., Padoan & Nordlund 1999). 

2.2. Generation of Cloudy ISM through Thermal Instability 

From the above-mentioned understanding of interstellar 
clouds, we employ a cloudy ISM formed by the thermal in- 
stability as an ambient medium of a SNR. We solve the three- 
dimensional ideal MHD equations with interstellar cooling, 
heating, and thermal conduction: 

|w.(pvl = 0, 

dpv -L ^ ^ ^ B^B^ ^ 
^+V(/7+ — ) = 0, 

Ot 8 TT 4 TT 

Oe B^ B V* 

— +V-{(e+p+—)v-—B} = V-KVT-pCin,T), 

at O TT 4 TT 

dB ^ ^ ^ 
-- = Vx(vxB), 
at 

p pv^ B^ 



where k, is the thermal conductivity and C{n^T) is the net 
cooling rate per unit mass. We employ a net cooling func- 
tion that is obtained by fitting various line-emission coolings 
(Ly-a, CII 158/im, 01 63 /im, etc.) and photoelectric heating 
by dust (Koyama & Inutsuka 2002), which can adequately de- 
scribe the effects of cooling and heating roughly in the tem- 
perature range 10K<r<10'^K. We impose the ideal gas 
equation of state and the adiabatic index 7 = 5/3. Since the 
thermal conduction determines the most unstable scale of the 
thermal instability (Field 1965), it should be taken into ac- 
count during cloud formation. Because the medium is in a 
weakly ionized state, the isotropic thermal conductivity due 
to the neutral atomic collisions n = 2.5 x 10^ erg cm~^ s~^ 
is used (Parker 1953). The numerical technique used to 
solve the basic equations is a combination of a second-order 
Godunov-type finite volume scheme (Sano et. al. 1999) and 
a second-order consistent method of characteristics with con- 
strained transport algorithm (Clarke 1996). The former tech- 
nique enables us to sharply capture a shock wave (van Leer 




Fig. 2. — Number density volume rendering of the resulting cloudy 
medium as a consequence of the thermal instability after 3.0 Myr of evo- 
lution (a few cooling times). The number density map in the _y = 0.0 pc plane 
is overplotted. Regions in green and blue indicate the density ^ ~ 10 cm~^ 
and n >30 cm~^, respectively, and the region in red shows the diffuse in- 
tercloud gas with n < 1 cm~^ . Magnetic field lines are represented as gray 
lines. 

1979), and the later allows us to integrate the induction equa- 
tion without breaking the divergence free condition of mag- 
netic field (Evans & Hawley 1988) and to stably follow the 
strong magnetic field amplification due to turbulence (Clarke 
1996). 

We use a cubic numerical domain whose side lengths are 
2 pc with the resolution of Ax = 2 pc/1024 = 1.95 x 10"^ pc 
in which periodic boundary conditions are imposed. As the 
initial conditions, we choose a thermally unstable equilibrium 
state with random density fluctuations whose mean density 
and thermal pressure are ^ = 2.0 cm~^ and ;?/^b = 2887 K 
cm~^, respectively. The uniform magnetic field oriented -\-y 
direction is imposed whose strength is 5.0 /iG, which is be- 
lieved to be the average strength in the ISM (Beck 2000). In- 
oue & Inutsuka (2008, 2009) showed that such a thermally 
unstable medium can be ubiquitously expected as an initial 
condition of interstellar clouds. 

Fig. [2] shows the number density volume rendering of the 
resulting cloudy medium as a consequence of the thermal in- 
stability after 3.0 Myr of evolution (a few cooling times). The 
number density map in the y = 0.0 pc plane is overplotted. 
Regions in green and blue indicate the density ^ ~ 10 cm~^ 
and ^ ^ 30 cm~^, respectively, and the region in red shows 
the diffuse intercloud gas with n ^ 1 cm~^. Magnetic field 
lines are represented as gray lines. The condensations driven 
by the thermal instability to form clouds arise along the mag- 
netic field lines, since the motion perpendicular to the field 
is easily stopped due to the magnetic pressure. This results 
in a formation of sheet-like clouds whose thicknesses are es- 
sentially determined by the most unstable scale of the thermal 
instability 1 pc) times the compression ratio of the con- 
densation (~ 0.1) that gives ~ 0.1 pc, and whose length is 
roughly given by the most unstable scale ~ 1 pc. Note that 
the intercloud gas is also thermally unstable due to runaway 
heating that evolves toward the diffuse gas phase. The typi- 
cal densities and temperatures of the clouds are nc ^ 40 cm~^ 
and Tc 100 K, and those of the diffuse intercloud gas are 
nd — l cm~^ and ^ 5,000 K. A more detailed description 
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of the evolution of the thermal instability can be found, e.g., 
in Inoue et al. (2007) and Inoue & Inutsuka (2008, 2009). 

In this simulation, the formed clouds correspond to HI 
clouds and their volume filling factor (where ^ > 10 cm~^) is 
1.9%. In §3 and §4, we discuss how cloud density and cloud 
filling factor affect the dynamics of SNR. We stress that even 
in the formation of denser molecular cloud, dense clumps are 
generated by the thermal instability and the clumps are em- 
bedded in the diffuse gas (see, e.g., Hennebelle et al. 2008; 
Banerjee et al. 2009). 



2.3. Induction of Shock Wave 

To study the formation of SNRs from the cloudy ISM, we 
induce a strong shock wave by setting a high-pressure hot 
gas with p\^/k^ = 1.0 x 10^ and % = 0.1 cm~^ at one of the 
boundaries of the simulation domain. We examine two cases 
of shock induction from the x = surface and the j = sur- 
face, which correspond to the perpendicular shock and the 
parallel shock, respectively. The resulting average propaga- 
tion speeds of the induced shock waves are 2429 km s~^ and 
2330 km s~^ for the perpendicular and parallel shock cases, 
respectively, which correspond to the supernova blast wave 
shocks with an age of ~ 10^ yr. For the perpendicular (paral- 
lel) shock case, we set the periodic boundary condition at the 
y {x) and z boundaries and the free boundary condition at the 
X = 2 (y = 2) pc boundary. This numerical setting is very simi- 
lar to the simulation performed in Inoue et al. (2010) in which 
only the perpendicular shock with a much smaller shock ve- 
locity (Vsh ~ 500 km s~^) is presented. 

In contrast to the stage of setting up the initial condition, 
we omit the effect of cooling in the simulation of shock prop- 
agation, since the cooling timescale in shocked diffuse gas is 
larger than the timescale of the shock crossing time. As for 
the cooling in the shocked cloud, it can also be also neglected 
in our choice of initial medium. According to the fitting line 
emission cooling rate given by Gaetz et al. (1987), the cooling 
timescale for shocked cloud is estimated to be 

^cooi:^3xlO ^3QQQj^^^-iJ U00cm-3j 
V 1 cm 3 / 

where Vsh is the shock speed in the diffuse medium and 
we have used the fact that the shock velocity in clouds is 
stalled by a factor of y/nj/nc (see, §3.1 below). This cool- 
ing timescale is longer than the age of typical young SNR 
< 1,000 yr. However, as we discuss in §4.3 and §4.4, the 
cloud density that seems to be interacting with the shock in 
RXJ1713.7-3946 may be much larger (typically nc > 1,000 
cm~^). In that case, the timescale of cooling becomes much 
smaller than the age of the SNR. As we have shown in Ap- 
pendix A, the dynamics outside the shocked cloud marginally 
depends on the effect of cooling in the shocked cloud. Thus, 
even if the effect of cooling in shocked cloud is very effec- 
tive (i.e., 7 = 1 in the cloud), we can expect the similar re- 
sults outside the very dense regions where generation of tur- 
bulence and magnetic field amplification take place (see §3.1 
and §3.2). We also omit the effect of thermal conduction, 
since the shocked gas becomes a fully ionized gas in which 
the gyro radius for thermal protons is much smaller than the 



resolution of our simulation: 




Note that, since the basic MHD equations are integrated 
in conservative fashion, our code can accurately handle high 
Mach number shocks. Since we perform this shock propa- 
gation simulations in the upstream-rest-flame, the kinetic and 
magnetic energy does not dominate the thermal energy every- 
where in the computational domain that enable us to perform 
the simulations of a very high Mach number shock propaga- 
tion even in the use of the total energy conservative code. 

3. RESULTS 
3.1. Generation of Turbulent SNR 

Propagation of the shock wave compresses and piles up the 
cloudy medium and forms a turbulent shocked slab. Figs. [3] 
and |4] show density and magnetic field strength structures of 
the perpendicular and the parallel shock cases at ^ = 750 yr af- 
ter the shock injection, respectively. When the shock wave 
hits a cloud, a transmitted shock wave propagates into the 
cloud whose propagation speed is slower than the shock ve- 
locity in the diffuse intercloud gas by a factor of ^/r^Jn'c (see 
Appendix A for details), where ni and n^ are the number den- 
sities of intercloud gas and cloud, respectively. The stalling 
of the shock wave can be easily understood as follows: In the 
case of a blast wave, the pressure behind the shock becomes 
roughly isobaric, except with fluctuations of an order of mag- 
nitude (see §3.4 below), and the shock velocity in the pre- 
shock rest frame is essentially determined by the post-shock 
sound speed. Thus, the shock propagation speed becomes in- 
versely proportional to the square root of pre- shock density. 
As a result of the stall, the shock front is deformed when it 
interacts with dense clumps. In the panels (b) of Figs. |3]and 
|4] the deformations of the shock waves are clearly observed. 
The deformation of the shock front due to a pre- shock den- 
sity fluctuation is known as the Richtmyer-Meshkov instabil- 
ity (see, Brouillette 2002, Nishihara et al. 2010 for reviews). 
The deformed shock wave leaves vorticity in the post-shock 
flow, even if the pre-shock is static (Crocco's theorem; see 
also Kida & Orszag 1990, Zabusky 1999) that drives turbu- 
lence. 

In addition to this curved shock effect, the baroclinic ef- 
fect due to the misalignment of pressure and density gradients 
also generates vorticity. In the simulations, the baroclinic ef- 
fect is significant at the transition layer between the cloud and 
diffuse gas where density gradient is very large. The vortic- 
ity generated by the baroclinic effect is the strong velocity 
shear flows between the shocked cloud and the diffuse gas. 
In other words, the velocity difference in the shocked cloud 
and shocked diffuse gas generates the velocity shear, since 
the shock speed is stalled in the cloud while the shock and 
post- shock gas around the cloud continue to propagate with a 
faster velocity. 

The velocity dispersion of the shocked gas for the perpen- 
dicular shock case is 705 km s~^ and that for the parallel shock 
case is 738 km s~^ The corresponding Mach numbers are, 
respectively, {M^Y^^=0.6\ and 0.62. Note that these disper- 
sions mainly represent the velocity dispersions of diffuse gas. 
Here we have defined the shocked gas as the gas with ^ > 0.5 
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Fig. 3. — Result of the perpendicular shock case at f = 750 yr after the shock injection. Panel (a): number density volume rendering. Regions in green and 
blue indicate the density /i ~ 10 cm~^ and n>30 cm~^, respectively, and the regions in warm colors show the shocked diffuse gas with n <4 cm~^. Magnetic 
field lines are represented as gray lines. Panel (b): two-dimensional number density slice at z = 0.0 pc. Panel (c): slice of number density atx = 1.5 pc. Panel (d): 
volume rendering of magnetic field strength. Regions in blue, green and red indicate the regions with 5 < 100 /iG, 5 > 100 /iG, and B > 500 /iG, respectively. 
Panel (e): slice of magnetic field strength at z = 0.0 pc. Panel (f): slice of magnetic field strength at x = 1.5 pc. 



cm~^ and /^/^b > 10^ K cm~^. The former condition excludes 
the thin hot gas injected from the boundary and the later con- 
dition excludes the pre- shock gas. In what follows, we use the 
above two conditions to select the shocked gas. 

As we mentioned in the previous section, most of the vol- 
ume in the initial medium is filled by diffuse gas. Thus, the 
global shock speed is essentially determined by the shock 
propagation speed in the diffuse gas. This indicates that even 
if a SNR is interacting with interstellar clouds, the global ex- 
pansion rate of SNR is determined by the shock propagation 
speed in the diffuse gas. 



3.2. Magnetic Field Amplification 

In the post-shock region, the magnetic field is strongly am- 
plified far beyond the shock compression value, which is very 
similar to the previous simulations (Inoue et al. 2009, 2010). 
The top panel in Fig. [5] shows the plots of the evolutions of 
maximum magnetic field strength and average field strength 
(1^1) in the shocked slab. Red and blue lines are results of 
the perpendicular and the parallel shock cases, respectively. 
In both the perpendicular and parallel shock cases, the max- 
imum magnetic field strength reaches on the order of 1 mG, 
which is 200 times larger than the initial strength. We also 
plot the evolution of plasma (3 = ^7Tp/B'^ where the magnetic 
field strength is maximum in the bottom panel of Fig. [5] It is 
clear that the saturation of the magnetic field amplification is 
determined by the condition of /3 ~ 1 . 

One can understand the mechanism of this strong amplifi- 
cation using the following equation that is obtained from the 



equation of continuity and the induction equation: 




(4) 



This equation indicates that the magnetic field is amplified, 
if the velocity has a shear along the field line. We can also 
express the amplification as Faraday's law of induction, since 
the velocity shear along the magnetic field line is equivalent 
to a rotation of the electric field. Such a situation is realized, 
in particular, at the transition layer between the cloud and dif- 
fuse gas due to the baroclinic effect. From equation (|4]), we 
can estimate the timescale of the magnetic field amplification 
around the clouds as ^grow ^ ^tr/Av, where /tr is the thickness 
of the transition layer between the cloud and diffuse gas and 
Av is a velocity difference of the post-shock gas flows be- 
tween the cloud and the diffuse gas (or a difference of the 
shock propagation speed). Since the shock speed in the cloud 
is much smaller than that in the diffuse gas, the velocity dif- 
ference Av can be estimated by the shock speed in the diffuse 
gas 2000 km s~^). The thickness of the transition layer 
is given by the so-called "Field length" /tr ~ /f = ^/hiT/pC, 
which expresses the thermal balance between structure forma- 
tion by cooling and diffusion due to conduction (Field 1965, 
Begelman & McKee 1990), where Hi is the thermal conductiv- 
ity and C(n, T) is the cooling rate per unit mass. In a typical 
ISM, /tr takes a value of 0.05 pc (Inoue et al. 2006). The above 
values give the growth timescale of ^grow ~ 20 yr. The results 
of simulations show that the growth times of the magnetic 
field at which the maximum magnetic field strength achieve 
/comp ^ini e (onc ^-folding time of the amplification after the 
shock passage, where /comp = 4 for the perpendicular shock 
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Fig. 4. — Same as Fig.|3]but for the parallel shock case. 
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Fig. 5. — Top panel shows the evolution of maximum magnetic field 
strength and average field strength in the shocked gas. Bottom panel shows 
evolution of plasma (3 where the magnetic field strength is maximum. Red 
and blue lines are the results of the perpendicular and the parallel shock cases, 
respectively. 

case and /comp = 1 for the parallel shock case) are 39 yr for 
the perpendicular case and 17 yr for the parallel shock case, 
which are roughly in agreement with the above estimation. 
Note that this strong amplification can be calculated correctly 



only if the scale of the transition layer ~ 0.05 pc is well re- 
solved. Our simulation with Ax c^2x 10~^ pc satisfactorily 
expresses the transition layer with more than 25 grid cells. 
Also note that the strength of the vorticity generated behind 
the shock in the surface region of a cloud is determined by 
the upstream Field length, because the baroclinic vorticity is 
caused by the interaction of upstream density gradient and the 
shock. Thus, the Field length mentioned above indicates the 
Field length in the upstream medium. 

So far, we have only discussed the magnetic field amplifica- 
tion at the transition layer of the clouds and diffuse gas. How- 
ever, amplifications also arise in the diffuse gas, because the 
average magnetic field strength {\B\) shown in Fig. [5] mainly 
reflects the magnetic field strength in the diffuse gas. In fact, 
in Fig. |2]and Fig. [3] we can confirm the amplification in re- 
gions other than the thin transition layers. As we have pointed 
out in |3.1[ turbulent eddies in the diffuse gas are induced by 
the effect of curved shock wave. The induced eddies amplify 
the magnetic field through the stretching of magnetic field 
lines, or, in other words, through the action of the turbulent 
(or small-scale) dynamo effect (Brandenburg & Subramanian 
2005, Cho & Vishniac 2000, Cho et al. 2009). Note that such 
an amplification as the result of interaction between a shock 
and density fluctuation has been recognized by many authors 
(Giacalone & Jokipii 2007, Inoue et al. 2009, 2010, 2011, 
Mizuno et al. 2010). Thus, as a consequence of the shock- 
cloud interaction, we can observe the two types of magnetic 
field amplification in response to the two mechanisms of vor- 
ticity generation. 

Prior to the recognition of the magnetic field amplification 
by the interaction between a shock and density fluctuation, 
Balsara et al. (2001, 2004) found that the interaction be- 
tween the shock and interstellar turbulence induce magnetic 
field amplification. By measuring the field line stretching and 
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Fig. 6. — Power spectra of velocity (solid) and magnetic field (dashed). 
Red and blue lines represent the perpendicular and the parallel shock cases, 
respectively. Spectra are calculated using the data of the turbulent regions 
of X G [1.0, 1.5] nye [1.3, 1.8] HzE [0.0,0.5] (for the perpendicular shock 
case) andxG [0.3,0.8] njG [1.1,1.6] nzG [0.9,1.4] (for the parallel shock 
case) at r = 750 yr. The black line shows the Kolmogorov law k~^^^ . 

the Lyapunov exponents of the chaotic, heHcal flows created 
by the shock-turbulence interaction, Balsara & Kim (2005) 
showed that stretch, twist and fold (STF; Childress & Gilbert 
1995) of magnetic field line plays a crucial role in the mag- 
netic field amplification. On the other hand, in the interac- 
tion between a shock and density fluctuation, we have some- 
what strong evidence that shows the STF is not playing a 
main role in the amplification, because the result of the pre- 
vious 2D simulations with almost the same initial condition, 
in which the STF is definitively unable to work, shows very 
similar magnetic field amplification (Inoue, Yamazaki, & In- 
utsuka 2009, see also similar 2D simulations by Giacalone 
& Jokipii 2007 and IMizuno et al. 2010). This indicates that 
the magnetic energy in the results of our simulations is not 
mainly contained in the mean-field due to the STF, but con- 
tained in waves and/or locally stretched magnetic field lines 
generated by vortical and shear flows. Note that this result 
does not contradict Balsara' s conclusion, because the simula- 
tions by Balsara et al. (2001, 2004) involves realistic turbu- 
lence in upstream ISM, while the simulations of the density 
fluctuation and shock interaction do not. We expect that if our 
initial medium involved turbulence that creates large helicity 
density by the shock-turbulence interaction as shown in Bal- 
sara et al. (2001, 2004), we might also obtain the additional 
mean-field amplification by the STF. This speculation will be 
tested in our future study of the interaction between turbulent 
molecular cloud and a supernova blast wave. 

3.3. Spectra of Turbulence 

Fig. |6] shows the one-dimensional power spectra of the 
velocity (Pvik)', solid) and magnetic field (PB(k); dashed), 
where the power spectra are defined as J Pydk = J d^x and 
J PBdk = J d^x. Red and blue lines represent the perpen- 
dicular and the parallel shock cases, respectively. Spectra 
are calculated using the data of the turbulent regions of x G 
[1.0, 1.5] n J G [1.3, 1.8] n z G [0.0,0.5] for the perpendicular 
shock case and x e [0.3,0.8] D y e [1.1,1.6] n z e [0.9, 1.4] 
for the parallel shock case at t = 750 yr. These regions are 
selected so that more than 85% of the regions are filled by 
the shocked gas. The obtained spectra are very similar to 
those of the super- Alfvenic turbulence with isotropic, large- 
scale, divergence-free forcing (Cho & Vishniac 2000, Cho & 
Lazarian 2003, Cho et al. 2009). In large scales (/ = 27r/k> 
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Fig. 7. — Pressure slice at the x = 1.675 pc plane of the perpendicular 
shock case at r = 750 yr. 

0.03 pc), the velocity field shows the Kolmogorov spectrum 
(X k~^^^, and the spectrum of the magnetic field is nearly flat. 
In small scales, the power spectra go down likely due to nu- 
merical dissipation. 

Theoretically, below the scale where the velocity disper- 
sion becomes comparable to the mean field Alfven velocity 
due to the cascade of the turbulence, the cascade mechanism 
can be changed from the Kolmogorov-like one to the IMHD 
critical balance cascade (Goldreich & Sridhar 1995). Since 
the velocity dispersion of Kolmogorov turbulence depends on 
the scale as /^/^, if we use the driving scale of turbulence 
L ~ 0.2 pc from Fig. [6j the velocity dispersion at the driv- 
ing scale {v)l ~ 700 km/s, and the mean field Alfven ve- 
locity (va) - 100 km s"^ (1(^)1/100 /iG)((n)/4cm-3)i/2, the 
transition scale of the cascade mechanism can be evaluated as 

/cas-Z.(vA)V(v)i ^5x10-4 pc. 

In the scales larger than /cas, we can expect large-amplitude 
isotropic magnetic field fluctuations, because the magnetic 
field is passive with respect to the turbulent velocity field. 
Thus, the resonant scattering of particles by Alfven waves 
would be very effective and we can expect the B ohm-limit 
diffusion (SB^/B^ 1), if the resonant scale of particles /j-es = 
lirE/eB 1 x 10'^ pc (E / lOOTcY) (B / 100 i^G)'^ is larger 
than /cas- In the scales smaller than ^as, according to Yan & 
Lazarian (2002), the resonant scattering of particles by Alfven 
waves would be ineffective due to the anisotropy of MUD tur- 
bulence. However, this does not indicate the inefficiency of 
the scattering by the IMHD turbulence for the particles whose 
resonant scale is smaller than /cas- Beresnyak et al. (2010) re- 
cently studied the transport of test particles in a compressible 
MHD turbulence. They found that the particles are effectively 
scattered by non-resonant mirror reflections due to the large- 
scale magnetic field fluctuations. Thus, the IMHD turbulence 
induced in the simulations would be important for the particle 
scattering even for the low-energy particles whose the reso- 
nant scale is smaller than /cas- 

3.4. Transmitted and Reflected Shock Waves 

When the strong shock wave hits a cloud, a transmitted 
shock wave penetrates into the cloud and a reflected shock 
wave propagates back into the shocked slab. Since the re- 
flected shocks are formed at the surface of the clouds, many 
reflected shock waves are generated in the SNR. For instance, 
in Fig. [7] we show the pressure slice of the perpendicular 
shock case in the y-z plane at x = 1.675 pc. We can see a 
number of discontinuous pressure jumps due to the reflected 



8 



T. INOUE ET. AL. 



Mach number : M 

^ 10° IQi 102 10^ 

0.01 I — — — 




pressure jump : Ap 

Fig. 8. — Probability distribution function of pressure jump Ap. Red and 
blue lines are the results of the perpendicular and the parallel shock cases, 
respectively. Upper horizontal axis denote the Mach number of shock waves 
converted by using eq. |5j. 

shocks in the SNR. In the following, we also describe the re- 
flected shocks as the secondary shocks to distinguish them 
from the primary forward shock wave. 

In Appendix A, we analytically evaluate the Mach number 
of the reflected shock using a formulation given by Miesch 
& Zweibel (1994). We find that the Mach number of the re- 
flected shock wave is a function o f the density ratio of the 
intercloud gas and cloud (see eq. |A3| ). In the case of our 
simulation, the ratio niju^^ 40 gives the Mach number of the 
secondary shocks as M 1.8. We also show in Appendix A 
that the Mach numbers of the shock wave in diffuse gas and 
the transmitted shock wave in clouds are almost identical. 

To evaluate the strength of the transmitted and reflected 
shocks in the simulations, we plot the probability distribution 
function (PDF) of the pressure jumps in Fig. [8] The PDF 
is calculated using the following procedure: We define the 
sphere whose radius is 20 times the numerical grid size and 
whose center is at each cell in the numerical domain. In each 
sphere, we calculate the minimum pressure /^min? the maxi- 
mum pressure /^max, and the maximum pressure gradient, ex- 
cept for the injected hot plasma with n < 0.5 cm~^. If the 
maximum pressure gradient is larger than the critical pres- 
sure gradient pcr = /7max/(5 Ax), i.e., the pressure fluctuation 
in the sphere is caused within the narrow region, we regard 
a shock wave as being in the sphere whose pressure jump is 
Ap = PmsLx/Pmin- Here the factor 5 in the expression of the 
critical pressure gradient is chosen from the fact that the high- 
resolution shock-capturing scheme used in this study requires 

5 grid points at most to express a shock wave. The bimodal 
distribution of the pressure jump in Fig. [8]shows the existence 
of the primary shock and secondary shocks. According to the 
Rankine-Hugoniot relation, the Mach number of a shock for 
the gas of 7 = 5/3 is related to the pressure jump as (Landau 

6 Lifshitz 1959): 

M(Ap)=^^ + ^Ap. (5) 

Substituting the peak of the PDF A/? = 3.85 that corresponds 
to the secondary shocks, we obtain = 1.81, which agrees 
well with the above-mentioned analytical evaluation M^l.^ 
(see. Appendix A for more detail). In the upper horizontal 
axis of Fig. [8] we also exhibit the Mach number of shock 
waves converted by using eq. ([5]). 

The peaks of the PDF with higher pressure jumps are com- 
posed of the pressure jumps due to the primary shock prop- 



agating in the diffuse gas and the transmitted shocks in the 
clouds. This also agrees well with the analytical evaluation 
that the Mach numbers of both the shock in the diffuse gas 
and the transmitted shock are almost identical, i.e., the pri- 
mary shock propagates so that the Mach number remains un- 
changed even after the encounter with the cloud. 

Note that even in the case of a different initial medium, e.g., 
a medium with larger clump density and larger clump filling 
factor, the Mach number of the secondary shocks would be 
affected marginally. This is because the Mach number of the 
reflected shock depends only weakly on the density ratio of 
the intercloud gas and cloud, and the Mach number converges 
to = >/5 in the limit of a large density ratio (see Appendix 
A). 

4. DISCUSSION: OBSERVATIONAL FEATURES 

4.1. Short-Time X-ray Variability 

From the observations of RX J17 13.7-3946 using the 
Chandra and Suzaku X-ray space telescopes, Uchiyama et 
al. (2007) discovered that the synchrotron emissions from X- 
ray hot spots whose spatial scale is ~ 0.05 pc show timescale 
variability of a few years. They concluded that the short- time 
variability would be a consequence of magnetic field ampli- 
fication up to ImG. This is because the timescale of the syn- 
chrotron cooling that causes the decay of X-ray luminosity, 
and the acceleration timescale of the B ohm-limit DSA that 
causes brightening, are respectively given by (Uchiyama et 
al. 2007, Malkov & Drury 2001) 

^--^Ht^J (t^) (3,oooLs-0 

(7) 

where r] = dB^/B^ is the degree of magnetic field fluctuations 
that characterizes the efficiency of the acceleration and e = 
1.6(5/lmG)(£'/10TeV)^ keV is the energy of synchrotron 
photons emitted from accelerated electrons with energy E. 

The results of our simulations can explain such a short-time 
variability of the synchrotron X-rays. As shown in §3.2, the 
strong magnetic field amplification easily makes regions with 
B ^ I mG, and their spatial scale ^ 0.05 pc also agrees quite 
well with the observed regions. Since our results show the 
regions of 5 ~ 1 mG are downstream of the primary shock 
wave, an additional accelerator other than the primary shock 
is necessary to produce brightening of the X-rays. The sec- 
ondary shocks generated by the shock-cloud interaction can 
be such accelerators. Although the Mach number of the sec- 
ondary shocks is much smaller than that of the primary shock, 
their propagation speed is comparable to the primary shock's, 
because they are downstream of the primary shock where the 
sound speed is comparable to the primary shock speed in a 
diffuse medium. Thus, Vsh ^ 3000 km s~^ in eq. Q is avail- 
able even when the accelerator is the secondary shock. Owing 
to the small Mach number of the secondary shocks, the injec- 
tion rate of the electron acceleration may be much smaller 
than the case of the primary shock. Even so, the secondary 
shocks can accelerate electrons, because the secondary shocks 
are in the downstream region of the primary shock where the 
relativistic electrons accelerated by and advected from the pri- 
mary shock are available as seed particles for further acceler- 
ation. Note that, in order to adopt the estimation of eq. ([7]), 
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both upstream and downstream magnetic field should be am- 
plified and such high-field region should be larger than the 
acceleration sites, i.e., the region in which the particles being 
accelerated exist. Again, since the secondary shocks propa- 
gate in the shocked region, the magnetic field is already am- 
plified upstream of the secondary shocks and electrons with 
an energy of 10 TeV may go upstream and reach the distance 
of rjlgc/vsh — 5 X lOr^T] pc from the shock front (assuming 
B mG), which is smaller than the scale of the milli-Gauss 
regions ~ 0.05 pc. Therefore, the acceleration timescale of a 
few years is basically possible at the secondary shocks. Note 
that the spectrum of particles reaccelerated by the secondary 
shock would be steeper than that accelerated by the primary 
shock because of the small Mach number of the secondary 
shocks. This indicates that the synchrotron emission of elec- 
trons reaccelerated by the secondary shocks can dominate the 
emission of electrons accelerated by the primary shock only 
in the regions with high magnetic field strength (i.e., the elec- 
trons accelerated by the primary shock are already cooled)^ . 

If we only take into account the resonant scattering of par- 
ticles by turbulent MHD waves, the B ohm-limit accelera- 
tion {j] = 1) in eq. ^ may not be achieved because the 
resonant scale for relativistic electrons I^q^^IitE / eB ^1 x 
10~^ (£^/10TeV)(5/lmG)~^ pc is somewhat smaller than the 
transition scale of the cascade mechanism (/cas ~ 5 x 10""^) 
above which the Bohm-limit diffusion can be expected (see 
§3.3). However, as discussed in §3.3, we can additionally 
expect the non-resonant mirror scatterings by the large-scale 
magnetic fluctuations that would make the Bohm-limit accel- 
eration possible in the turbulent shell. Furthermore, accord- 
ing to Ohira et al. (2009), the effect of a cold proton beam 
created by a charge-exchange process in shocked neutral gas 
causes large-amplitude magnetic field fluctuations through the 
growth of the Weibel instability or the resonant instability that 
also enhances the acceleration efficiency around the shocked 
clouds. It is worth noting again that our scenario explains 
not only the timescale of X-ray variability but also its typical 
spatial scale of 0.05 pc, since the strong magnetic field am- 
plifications to milli-Gauss take place at the transition layers 
between the clouds and diffuse gas that is essentially deter- 
mined by the Field length ~ 0.05 pc (Inoue et al. 2006). 

The evidence of magnetic field amplification in RX 
J 17 13.7-3946 has also been obtained from the thickness of 
X-ray filaments that typically show 5 ~ 100 fiG (Hiraga et 
al. 2005, Ballet 2006). Since the X-ray filaments have ap- 
parently coherent features on parsec scales, it seems to be the 
average field strength. The average field strengths in our sim- 
ulations are tens of micro-Gauss (see Fig. [5]), which possibly 
explains the X-ray filaments. However, the average strength 
in our simulations becomes noticeable only in the downstream 
region far from the primary shock front suggesting that other 
possible amplification mechanisms such as the cosmic-ray 
streaming instability (Lucek & Bell 2000) and/or dynamo ef- 
fect driven by the interaction between cosmic-ray precursor 
and density fluctuation (Beresnyak et al. 2009) might be im- 
portant in the vicinity of the primary shock front. We also 
stress that as pointed out by Pohl et al. (2005), the evaluation 

^ In the case of a middle-aged SNR, since the acceleration by the primary 
shock can be ineffective above the break energy of ~ 10 GeV due to the 
damping of MHD waves by ion-neutral collisions, the particles accelerated 
by the secondary shocks can dominate the spectrum above the break energy 
(see, Inoue, Yamazaki, & Inutsuka 2010 for detail). This is a possible origin 
of the broken power-law spectrum of particles commonly observed in middle- 
aged SNRs interacting with molecular clouds. 



of the magnetic field strength in the X-ray filaments based 
on the synchrotron cooling timescale might be questionable, 
i.e., the thickness of the filaments might be determined by the 
spatial scale of the strong magnetic-field region instead of the 
cooling length of electrons in a uniformly magnetized region. 

4.2. Global Morphology of SNR 

One may wonder that despite its associating molecular 
clouds, especially in the northeast region, the morphology of 
RX J17 13.7-3946 is roughly spherical and the expansion rate 
is fast as if it is evolving in the diffuse environment. This is 
not surprising if we take into account the clumpy nature of 
the molecular clouds: As we showed in §3.1, the mean prop- 
agation speed of the primary shock is determined by that in 
the diffuse gas. The speed of the primary shock is stalled lo- 
cally where it hits a clump of the clouds. However, once the 
shock passes the clump or once the shock entirely surrounds 
the clump and is topologically separated from the clump, the 
shock is accelerated and the deformation is dissolved, because 
the fast shock is stable with respect to the deformation. Thus, 
the local configuration of the shock wave can be deformed 
only where the shock is interacting with local clumps, while 
the global configuration would be affected only if the density 
of the diffuse surrounding of the clouds has global variations. 

4.3. Suppression of Thermal X-ray Line Emission 

It has been shown that RX J17 13.7-3946 is characterized 
by the lack of thermal X-ray emission (Slane et al. 1999, 
Cassam-Chenai et al. 2004, Takahashi et al. 2008). Ac- 
cording to Takahashi et al. (2008), the upper limit of gas 
density to ensure the thermal X-ray-less feature is 0.01 to 2 
cm~^ depending on the assumed electron temperature. The 
upperlimit density {n ^ 0.01-2 cm~^) is smaller than the mean 
density of the diffuse ISM (^ ~ 1 cm~^), if the electron tem- 
perature is larger than 0.05 keV. In addition, if the temperature 
of shocked clouds is larger than ~ 1 keV, they would emit an 
unacceptable amount of thermal X-ray lines because of their 
high density. However, if we take into account the effect of 
stellar wind from the progenitor of a supernova, the density of 
the diffuse gas can be smaller than the threshold density for 
thermal X-ray line emission. Furthermore, the thermal X-ray 
line emission from the shocked clouds are substantially sup- 
pressed because of the stall of the transmitted shock. These 
two effects on the thermal X-ray line emission were pointed 
out in Zirakashvili & Aharonian (2010). In the following, we 
discuss these effects more quantitatively based on the results 
of our simulations. 

The requirement for the density of the diffuse gas can be 
achieved if the progenitor of RX J17 13.7-3946 is a massive 
star as is widely believed (Slane et al. 1999). This is be- 
cause the stellar wind from the massive star would sweep up 
pre-existing intercloud gas rarefying the intercloud gas signif- 
icantly, while dense clumps are not swept off owing to their 
high density (e.g., Gritschneder et al. 2009). The situation is 
illustrated schematically in Fig. [9] According to Weaver et 
al. (1977), who studied the expansion of a bubble formed by 
stellar wind from 0-type stars, the resulting gas density in the 
wind bubble is ^ ~ 0.01 cm~^ (see, e.g.. Fig. 3 of Weaver et al. 
1977). Note that the density in the wind bubble is not deter- 
mined by the density of wind gas but by the evaporation of the 
wind shell into the cavity. The radius of stellar wind bubble 
is described using the mechanical luminosity of the wind 
Lw, density of interstellar gas no, and lifetime of the wind 
%e as K = 27 pc (Lw/lO^^ ,-i)i/5 ^^^/^ cm-3)-i/5 ^t,^ji 
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Fig. 9. — Schematic view of wind bubble expanding in a cloudy ISM. 
Diffuse intercloud gas is swept by the stellar wind, while dense cloud cores 
and clumps can survive in the wind. Density in the wind bubble is much 
smaller than the intercloud gas density that is determined by the evaporation 
of the wind shell by thermal conduction. 



Myr)^/^ (Castor et al. 1975). According to this expansion law, 
in order for the dense gas to stay within the cavity of the wind 
bubble, the density should be at least larger than 



^0 ^ 10 cm 



10^6 erg s-i ) V 10 pc 



1 Myr 



, (8) 



where we have adopted a distance of 1 kpc and thus the radius 
of RX J1713.7-3946 - 10 pc (Fukui et al. 2003). 

Recently, Sano et al. (2010) have shown by using the NAN- 
TEN telescope that the "peak C" of a CO molecular cloud 
core associated with the region in RX J17 13.7-3946 seems to 
be embedded in the SNR. Since the density of the molecular 
cloud core is approximately 10"^ cm~^, it is reasonable for such 
a dense object to stay in the SNR. Eq. ([5]) suggests that less 
dense molecular cloud cores or molecular clumps with den- 
sity on the order of 10^ cm~^ depending on and ^age would 
also be embedded in RX J17 13.7-3946, although these may 
not be observed by CO line emission surveys due to the dis- 
sociation of molecules by UV radiations from the progenitor 
massive star. We conclude that if we take into account the 
effect of the stellar wind from the massive progenitor, the dif- 
fuse intercloud gas density becomes on the order of ^ ~ 0.01 
cm~^, which does not conflict the lack of the thermal X-ray 
line emission, while dense molecular clumps/cores can be left 
in the wind bubble. 

The remaining issue for the X-ray line emission from the 
shocked clouds is resolved easily as follows: The tempera- 
ture of protons in the shocked gas, which corresponds to the 
maximum temperature of electrons, is given by 



16"^ 



= 18 



Vsh 



3,000kms- 



keV, 



(9) 



where Vsh is the shock velocity that is supposed to be 3000 km 
s~^ in the diffuse gas (gas in the wind cavity with the density 
~ 0.01 cm~^). In the cloudy ISM, the shock is stalled when 
it hits a cloud. As we show in §3.1 and Appendix A in more 
detail, the velocity ratio of the shock wave in the diffuse gas 
and the cloud is proportional to the square root of their den- 
sity ratio: Vsh,d/vsh,c — (^c/^d)^^^- From this relation, we can 
estimate the proton temperature (corresponding to the upper 



bound of the electron temperature) of the shocked cloud as 
3 

16" 



= 2X 10" 



7 VO.Olcm-V 



VI 03 cm- 



3,000kms-i J VO.Olcm-^ 
1 

keV. 



(10) 



Therefore, even after the passage of the shock wave in the 
clouds, bright thermal X-ray line emission from the clouds is 
not expected. 

4.4. Spectrum ofHadronic Gamma Rays 

Recently, using one-dimensional model assuming a uni- 
form ISM, Ellison et al. (2010) claimed that if we reduce 
the ambient density to reconcile the absence of the thermal 
X-ray line emission from RX J17 13.7-3946, the hadronic 
gamma-ray emission becomes dim owing to the low target- 
gas density for tt^ creation. The reason is as follows: Ac- 
cording to Aharonian et al. (2006), the total gamma-ray en- 
ergy measured from 0.2 to 40 TeV in RX J1713.7-3946 is 
W 6 X 10"^^ (J/1 kpc)^(^tg/l cm"^)"^ erg, where J is a dis- 
tance and ^tg is a mean target gas density. Thus, supposing the 
low-density ISM, the efficiency of particle acceleration be- 
comes 100(^tg/0.06 cm"3)(^/10^^ erg)% indicating that the 
hadronic gamma-ray emission cannot be as bright as observed 
even if the acceleration is extremely efficient. 

However, in our shock-cloud interaction model, the 
hadronic emission from the clouds embedded in the SNR 
can be expected, because the high density shocked clouds do 
not emit thermal X-ray lines owing to the low-temperature as 
shown in eq. (TTO]). If we assume a typical density of clumps 
Hci ^10^ cm~^~and their volume filling factor / ~ 10"^, the 
effective mean target density can be rewritten as ^tg — ^ci/ 
and thus the efficiency becomes 6(^ci/10^ cm~3)(//10~3)%. 
Although precise evaluation of the filling factor / is beyond 
the scope of this paper, our model can reproduce the hadronic 
gamma-ray emission that is compatible with the canonical ac- 
celeration efficiency ~ 10%. 

In the case of a uniform ISM model, the spectral energy 
distribution of the hadronic gamma rays directly reflects that 
of the accelerated nuclei roughly above the critical energy of 
the TT^ creation (~ 0.1 GeV), i.e., the photon index of the 
hadronic gamma-ray emission is = 2 for the standard DSA 
scenario. However, in our shock-cloud interaction model, as 
we discuss in the following, the spectrum may deviate from 
this conventional spectrum. Because of the heavy shock stall 
in the clouds, particle acceleration at the transmitted shock 
wave would be inefficient and most of the particles could be 
accelerated at the primary shock in the diffuse gas. The high- 
energy nuclei accelerated at the primary shock interpenetrate 
diffusively to nearby shocked dense clumps. The amount of 
the neutral pions generated through the collisions of acceler- 
ated nuclei and matters in the clouds, which eventually decay 
and emit gamma rays (Issa & Wolfendale 1981, Naito & Taka- 
hara 1994, Aharonian et al. 1994), is proportional to the mass 
of the dense clumps illuminated by the high-energy nuclei. 
Since the penetration depth into the clouds would depend on 
the energy of accelerated particles, the hadronic gamma-ray 
spectrum can be deviated from that of the conventional model 
supposing uniform ISM (see also Zirakashvili & Aharonian 
2010). 
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In our simulations, the clouds formed by the thermal insta- 
bility have a sheet-like structure. Using the penetration depth 
/pd of the accelerated protons into the cloud, the mass of the 
cloud illuminated by the high-energy protons M /pd p, 
where p 40 cm~^ is the density of clouds and R is the scale 
of the cloud sheet ~ 1 pc. This indicates that the mass of the 
interpenetrated region of the cloud is proportional to the pen- 
etration depth as long as the penetration depth is smaller than 
the thickness of the cloud sheet 2^ 0.1 pc. However, as dis- 
cussed in §4.3, low-density cloud envelope would be wiped 
out by the stellar wind. In the case of RX J1713.7-3946, 
the clouds that remain in the SNR would be dense molecu- 
lar cloud cores/clumps as suggested in eq. ^ and observa- 
tionally in Sano et al. (2010). If the dense molecular cloud 
cores are gravitationally bound, collapsing objects, the den- 
sity structure becomes p{r) = po(r/ro)~^ (Larson 1969, Pen- 
sion 1969, Masunaga & Inutsuka 2000, Andre et al. 2000). 
This is indeed the case of the "peak C" of CO core embedded 
in RX J1713.7-3946 (Sano et al. 2010). In this case the mass 
illuminated by the high-energy protons can be written as 

M= / 47rr^po(r/ro)^dr = 47rporllpa (11) 

If the cores are gravitationally unbound, pressure-confined 
objects, the density would be approximately constant (Bon- 
ner 1956, Ebert 1955, Alves et al. 2001), and the illuminated 
mass can be written as 

M= [ 4^rV^r=^(37?2/pd-37?/p2, + /p^,) 

c^47rp7?2/pd for7?>/pd, (12) 

The above examples indicate M oc /pd as long as /pd <C 7?. The 
radius of the molecular cloud cores is typically 0.1-1 pc, and 
the scale of dense clumps in molecular clouds is usually larger 
than 0.1 pc because 0.1 pc is the minimum scale of the ther- 
mal instability (Field 1965, Koyama & Inutsuka 2004, Inoue 
et al. 2006) and the turbulent flows that form clumps by shock 
compressions is expected to be subsonic below 0.1 pc (Lar- 
son 1981, Heyer & Brunt 2004, Mac Low & Klessen 2004). 
Note that when the density of clouds are as high as 10^ cm~^, 
the cooling timescale of clouds becomes much smaller than 
the age of RXJ1713.7-3946 -1,000 years that substantially 
diminish the scale of shocked region of clouds due to cooling 
contraction. However, owing to the slowdown of shock speed, 
the penetration length of the shock wave into cloud in 1,000 
year is evaluated to be —0.01 pc (where we have assumed the 
density ratio between the clouds and diffuse gas to be 10^ as 
we discussed in §4.3) that is much smaller than the scale of 
the molecular cloud cores ^0.1 pc. Thus, we can still assume 
the minimum scale of clouds to be 0.1 pc, even if clouds are 
apparently embedded in SNR. 

Since the diffusion coefficient for the high-energy particles 
can be written as K,d = 4r]lgc/37T (Skilling 1975), where /g is 
the gyro radius of relativistic particles, the penetration depth 
due to the random walk of the particles can be written as 

= (lolv) (15^) [w^) 

where E is the particle energy, and ^age is the age of the SNR. 
For the magnetic field strength in the above expression, we 



have used the observed typical strength in the dense regions 
of molecular clouds with density — 10^ to 10^ cm~^ (Crutcher 
1999). The parameter t] = /6B^ has large ambiguity. As 
we have discussed in §3.3 and §4.3, the Bohm- limit diffusion 
can be expected not only around the shocked clouds but also 
in shocked clouds through the scattering by turbulent mag- 
netic field fluctuations (Beresnyak et al. 2011) and fluctua- 
tions generated by the effect of a cold ion beam (Ohira et al. 
2009). Also, the observations by Uchiyama et al. (2007) sug- 
gest that 7^ — 1 is realized at least around the clouds. Thus, 
if we assume 7^ ~ 1, the penetration depth for particles with 
£^ < 10 TeV can be smaller than the minimum scale of clouds 
~ 0.1 pc and the interpenetrated mass can be proportional to 
the square root of particle energy MiE) ex /pd(£') (xE^^^. 

In the conventional one-zone model that assumes the 
amount of target matter creating tt^ is independent of the en- 
ergy of accelerated protons, the spectral energy distribution 
of gamma rays becomes NiE) dE oc E~p dE above the critical 
energy for the tt^ creation and below the maximum energy of 
accelerated nuclei. Here p is the spectral index of the distri- 
bution of high-energy nuclei with p = 2 m the conventional 
DSA theory. While in the shock-cloud interaction model, the 
amount of matter creating tt^ depends on the square root of the 
energy of accelerated protons, so that the gamma-ray distri- 
bution becomes N^{E)dE ex MiE)E~P dE cx 1^^{E)E~p dE cx 
£-p+i/2 yields N^iE) cx E~^-^ for p = 2\ that is consis- 

tent with the recent observation of RX J1713.7-3946 (Abdo 
et al. 2011). This photon index of the hadronic gamma-ray 
emission (/?- 1/2) is the same as that of the inverse Compton 
emission (/?+ 1)/2 when p = 2. Thus, the spectra in the two 
scenarios are indistinguishable. Fortunately, as we discuss be- 
low, these two emissions can be distinguished if we focus on 
the spatial distribution of gamma rays. 

4.5. Spatial Inhomogeneity of Nonthermal Emissions 

In the previous section, we have shown that the emission 
mechanism of gamma rays cannot be distinguished from the 
gamma-ray observation alone. In the following, we discuss 
how we can clarify the emission mechanism. A significant 
difference between our shock-cloud interaction model and the 
conventional uniform ISM model is the spatial inhomogene- 
ity. The uniform ISM model predicts that emissions are spa- 
tially correlated from microwave to gamma ray irrespective of 
emission mechanism of gamma rays. On the other hand, the 
shock-cloud interaction model predicts the following charac- 
teristics (l)-(5). The schematic picture of the shock-cloud in- 
teraction model is given in Fig. [T0| 

(1) Synchrotron emission will be more powerful in the 
cloud-rich regions than the regions without clouds because 
of the turbulent amplification of the magnetic field as a con- 
sequence of the shock-cloud interaction. In other words, on 
several parsec scales, the X-ray emission will spatially corre- 
late with the CO distribution. Note that there is no necessity 
to always find CO emission near the X-ray bright regions, be- 
cause a considerable amount of CO molecules in clouds can 
be dissociated due to the UV radiation from the massive pro- 
genitor star and also because shock heating dissociates CO 
molecules in the shocked clouds. 

(2) In small scales on the order of sub parsecs, the local 
peaks of X-ray emission will show anti-correlation with the 
local peaks of CO emission, because the magnetic field ampli- 
fications arise most strongly around the clouds, in particular, 
at the transition layers between the clouds and diffuse gas. 
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Fig. 10. — Schematic picture of the shock- cloud interaction model. Pri- 
mary forward shock wave propagates through the cloudy wind bubble, where 
particle acceleration operates. Transmitted shock waves in clouds are stalled, 
which suppresses thermal X-ray line emission and particle acceleration in 
clouds. Shock-cloud interactions induce shock deformations and turbulent 
eddies. The turbulent dynamo effect amplifies the magnetic field that en- 
hances synchrotron emission. Secondary reflected shock waves are generated 
when the primary shock hits clouds that induce the short-time variability of 
synchrotron X rays where magnetic field strength is ~ 1 mG around shocked 
clouds. Hadronic gamma rays are emitted from dense clouds illuminated by 
accelerated protons whose photon index can be p- 1/2 = 1.5 for p = 2. 

(3) Since the magnetic field strength maximally grows to 
the level of 1 mG near clouds, the short- time variability of X- 
rays can be found in the X-ray bright regions especially in the 
vicinity of the clouds (see, §4.1 for detail). 

(4) The primary shock wave propagates with high velocity 
in the diffuse gas where synchrotron filaments can be formed 
as the SNRs in diffuse circumstances (e.g., Vink & Laming 
2003, Bamba et al. 2003, 2005) and leptonic gamma-ray 
emission would also be emitted, while in the clouds the trans- 
mitted shock waves are stalled where the particle acceleration 
is inefficient. 

(5) If the hadronic gamma rays emitted from clouds are 
more powerful than the leptonic emission from the primary 
shock in the diffuse gas, the distribution of gamma-ray emis- 
sion will show good spatial correlation with CO distribution. 
Note again that there is no necessity to always find CO bright 
regions at the gamma-ray bright regions because of the CO 
dissociation by UV radiation and shock heating. 

In the case of RX J1713.7-3946, the distribution of bright 
X-ray regions are globally well correlated with the distribu- 
tion of CO line emission, which is consistent with the feature 

(1) . In addition, recent observation by Sano et al. (2010) 
showed that the local peaks of X-rays are located around the 
local peaks of CO line emission, which supports the feature 

(2) of our scenario. Furthermore, the regions that show the 
short-time variability of X-rays discovered by Uchiyama et al. 
(2007) are located in the CO rich region that is also consistent 
with the feature (3). 

Although the spatial resolution of gamma rays in RX 
J17 13.7-3946 is not sufficient to be compared with the dis- 
tribution of CO line emission, these two seem correlated in 
the sense that the gamma-ray emissions are stronger in the 
north-west region where the CO emission is also strong that 



is consistent with the feature (5) of our model. Future gamma- 
ray observations with higher spatial resolution may clarify the 
correlation. However, in the southeast region, gamma rays are 
detected despite there is no CO emission. Since the flux of 
hadronic gamma-ray emission depends linearly on the mass 
of clouds, this may suggest two possible interpretations of 
the gamma-ray emission in the southeast region. One is that 
CO molecules are dissociated by UV radiation because of a 
smaller column density of clouds than other regions or CO 
molecules are dissociated by shock heating due to the smaller 
column density than other regions. In that case, HI emission 
may be found to compensate for the missing mass. Another 
possibility is that there are no clouds in the southeast region, 
and gamma-ray flux from that region is determined by the lep- 
tonic emission. Recent observation strongly support the for- 
mer possibility: Fukui et al. (2011) analyzed both HI and CO 
in RXJ1713 and have shown that the total ISM protons in both 
molecular and atomic phases exhibit a good spatial correlation 
with the distribution of TeV gamma-rays. It is noteworthy that 
the southeastern part of the gamma-ray shell shows an ISM 
counterpart only in the HI seen as self-absorption, while the 
rest of the shell are associated with both CO and HI. This cor- 
relation indicates that the ISM protons are a reasonable can- 
didate for target protons in the hadronic interaction, favoring 
that the hadronic interaction plays a role in the production of 
gamma-rays. 

5. SUMMARY 

We have examined the propagation of a strong shock wave 
(Vsh ~ 2500 km/s), which corresponds to a supernova blast 
wave shock with the age of ~ 10^ years, through a cloudy 
medium formed as a consequence of the thermal instability by 
using three-dimensional MHD simulations. We found that the 
shock-cloud interaction leads to deformation of a shock front 
and leaves vortices or turbulence behind the shock wave. The 
magnetic field behind the shock wave is amplified as a result 
of the turbulent dynamo action. The maximum magnetic field 
strength reaches up to 1 mG that is determined by the con- 
dition of plasma /3 ~ 1 . This is consistent with the previous 
simulations performed in limited two-dimensional geometry 
(Inoue et al. 2009). The scale of the region where B \ mG 
is determined by the thickness of the transition layer between 
the cloud and surrounding diffuse gas (~ 0.05 pc) at which 
the vortex is induced most strongly. The shock-cloud inter- 
actions generate many secondary shocks in the SNR at which 
particle acceleration can operate. The acceleration due to the 
secondary shocks in the region with B \ mG would be the 
origin of the short-time variability of X-rays discovered in the 
SNR RX J1713.7-3946 (Uchiyama et al. 2007). 

We also found that, since the medium formed as a con- 
sequence of the thermal instability is very clumpy, a shock 
wave propagating in the cloud is stalled heavily, while a shock 
in the diffuse gas is not. This gives the following important 
features of SNRs interacting with interstellar clouds: (1) The 
global morphology of the SNR interacting with clouds is not 
substantially affected by the clouds, since it is determined by 
the shock wave propagating in the intercloud, diffuse gas that 
fills the most volume. Thus, we should bear in mind that the 
shock velocity measured by its expansion rate is the velocity 
in the diffuse gas, and the shock velocity propagating in the 
cloud is much smaller. (2) The temperature of the shocked 
cloud is much smaller than the temperature in the post-shock 
diffuse gas (see eq.|[TQ|) that completely suppresses the ther- 
mal X-ray line emission from clouds. This could explain why 
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the thermal X-ray radiations from RX J 17 13.7-3946 are faint, 
despite the suggested interaction with molecular clouds. (3) 
Since the particle acceleration at the transmitted shock wave 
in clouds is inefficient due to the small shock velocity, the 
hadronic gamma rays from the clouds are emitted as a conse- 
quence of diffusion of high-energy protons accelerated at the 
shock wave in the diffuse gas. The penetration depths of the 
high-energy protons into the clouds depend on the square root 
of their energy that leads the photon index of hadronic gamma 
rays to be 1.5 (N(E)dE oc E-p^^^^ dE = E'^-^ for;? = 2). Thus, 
it is difficult to definitively distinguish the hadronic and lep- 
tonic gamma rays from the spectral study alone. We pro- 
pose that the detailed comparisons of the spatial distribution 
of gamma-ray. X-ray, and CO line emissions as discussed in 



§4.4 can be a conclusive method to reveal the origin of gamma 
rays from young SNRs. 
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APPENDIX 

A. STRENGTH OF TRANSMITTED AND REFLECTED SHOCKS 

When a SNR shock wave hits a cloud, a transmitted shock wave propagates into the cloud and a reflected shock wave propagates 
back into the shocked gas in the SNR. The situation is illustrated in Fig. |A1[ where the subscripts c, i, and s denote the values in 
the cloud, intercloud gas, and shocked intercloud gas, respectively. The symbols Vsh, Vr, and Vt represent the shock velocity in the 
diffuse gas, the transmitted shock velocity in the cloud, and reflected shock velocity measured in the rest-frame of the unshocked 
cloud (also the rest-flame of unshocked intercloud gas), respectively. An analytic treatment of this problem in plane parallel 
geometry is given by Miesch & Zweibel (1994). Owing to the facts that the primary shock is driven by thermal pressure, which 
acts isotropically, and the transmitted shock wave is stalled heavily due to the large cloud density, the primary shock tends to 
hits the cloud perpendicular to its surface. Thus, the following formula obtained by assuming one-dimensional geometry enables 
us to evaluate the typical strength of the shock waves, even though clouds have a complex structure. According to Miesch & 
Zweibel (1994), the Mach numbers of the transmitted and reflected shock waves are obtained by solving the following polynomial 
equations: 



27eM2 = ri^(27,>1?-7i + l)+7c-l, 
7i + l 

1 /2 

VTi/ 7c + l 7i + l 



(Al) 
(A2) 



where Mt = Vt/cc is the Mach number of the transmitted shock, Mr = (Vs - Vr)/cs is the Mach number of the reflected shock, 
Ms = Vs/cs is the Mach number of the shocked gas in SNR, S = pc/ps is the pressure ratio of the cloud and shocked gas, and 
e = ps/Pc is the density ratio of the shocked gas and cloud (again, the subscripts c, i, and s denote the values in cloud, intercloud 
gas, and shocked intercloud gas, respectively). 

By substituting eq. ( |A1[ ) into eq. ( |A2[ ), we can obtain a polynomial for Mr- In the case of the SNR interacting with cloud, S is 
a very small parameter. The series expansion of the polynomial for Mr with respect to S yields 



Mt 



7i + l 



MsMr- 



11 -e^M^ 
7c + l 



Mi 



1 



27: 



+ O((5) = 0. 



(A3) 



Since e is also a small parameter (but not very small compared to (5), Mr depends weakly on e. In the present case, 7i = 5/3, while 
7c can takes values from 1 (when the cloud density is sufficiently large to instantly cool down the shocked cloud) to 5/3 (when 
the cooling is inefficient). According to the Rankine-Hugoniot relation for a strong shock, the Mach number of the shocked gas 
in the rest frame of the pre-shock medium is given by M, = Vs/cs = {2/(7^ -7i)}^/^ (Landau & Lifshitz 1959). Thus, in the 
limit of e ^ 0, i.e., when the cloud is very dense, the physical solution of eq. dA3k gives Mr^VS^ 2.24 for 71 = 5/3, which 
corresponds to the upper-limit of Mr- This limit is independent of 7c. In our simulation (71 = 7^ = 5/3), the typical density ratio 
e ^ 0.1 gives Mr = 1.80. When 71 = 5/3 and 7c = 1, the solution is Mr = 1.75 indicating that the dynamics of the secondary 
shocks only marginally depend on the effect of cooling in the shocked cloud. 
From eq. (|A1|), the velocity of the transmitted shock wave into the cloud is approximately given by 



/ (7c + 1)7/ A 
7/ + 1 Pc 



Mr- 



(A4) 



On the other hand, from the Rankine-Higoniot relation for a strong shock, the shock velocity in the intercloud gas is Vsh = 
(7i + 1) Cs/{27i (7i - 1)}^/^. Thus, the ratio of the shock speeds in the cloud and intercloud is written as 



Vsh 



/ 27i(7c + l) Pi 
(7i + l)' Pc 




(A5) 
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Fig. A1. — Schematic diagram of the shock-cloud encounter in plane parallel geometry. The initial cloud and intercloud gas are in pressure equilibrium 
{pc = Pi). 

indicating that the propagation speed of the shock wave in cloud is slower than that in the shock in intercloud gas approximately 
by a factor of (pi/pc)^^^- Because the factor (pi/ Pc)^^^ is the ratio of the sound speeds in the cloud and diffuse gas, eq. ( A5 ) also 
indicates that thir Mach numbers are nearly equal. 
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